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NUMERICAL SCHEMES FOR KINETIC EQUATIONS IN THE 
DIFFUSION AND ANOMALOUS DIFFUSION LIMITS. PART I: THE 
CASE OF HEAVY-TAILED EQUILIBRIUM 

NICOLAS CROUSEILLES tt, HELENE HIVERT AND MOHAMMED LEMOU 

Abstract. In this work, we propose some numerical schemes for linear kinetic equations in the 
diffusion and anomalous diffusion limit. When the equilibrium distribution function is a Maxwellian 
distribution, it is well known that for an appropriate time scale, the small mean free path limit 
gives rise to a diffusion type equation. However, when a heavy-tailed distribution is considered, 
another time scale is required and the small mean free path limit leads to a fractional anomalous 
diffusion equation. Our aim is to develop numerical schemes for the original kinetic model which 
works for the different regimes, without being restricted by stability conditions of standard explicit 
time integrators. First, we propose some numerical schemes for the diffusion asymptotics; then, their 
extension to the anomalous diffusion limit is studied. In this case, it is crucial to capture the effect 
of the large velocities of the heavy-tailed equilibrium, so that some important transformations of the 
schemes derived for the diffusion asymptotics are needed. As a result, we obtain numerical schemes 
which enjoy the Asymptotic Preserving property in the anomalous diffusion limit, that is: they do 
not suffer from the restriction on the time step and they degenerate towards the fractional diffusion 
limit when the mean free path goes to zero. We also numerically investigate the uniform accuracy 
and construct a class of numerical schemes satisfying this property. Finally, the efficiency of the 
different numerical schemes is shown through numerical experiments. 

Key words. BGK equation. Diffusion limit, Anomalous diffusion equation, Asymptotic pre¬ 
serving scheme. 

AMS subject classifications. 35B25, 41A60, 65L04, 65M22. 


1. Introduction. The modeling and numerical simulation of particles systems 
is a very active field of research. Indeed, they provide the basis for applications 
in neutron transport, thermal radiation, medical imaging or rarefied gas dynamics. 
According to the physical context, particles systems can be described at different 
scales. When the mean free path of the particles (be. the crossed distance between two 
collisions) is large compared to typical macroscopic length, the system is described at 
a microscopic level by kinetic equations; kinetic equations consider the time evolution 
of a distribution function which gives the probability of a particle to be at a given 
state in the six dimensional phase space at a given time. Conversely, when the mean 
free path is small, a macroscopic description (such as diffusion or fluid equations) can 
be used. It makes evolve macroscopic quantities which depends only on time and on 
the three dimensional spatial variable. In some situations, this description can be 
sufficient and leads to faster numerical simulations. 

Mathematically, the passage from kinetic to macroscopic models is performed 
by asymptotic analysis. From a numerical point of view, considering a small mean 
free path, the kinetic equation then contains stiff terms which make the numerical 
simulations very expensive for stability reasons. In fact, a typical example is the 
presence of multiple spatial and temporal scales which intervene in different positions 
and at different times. These behaviors make the construction of efficient numerical 
methods a real challenge. 
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In this work, we are interested in the time evolution of the distribution function 
/ which depends on the time t > 0, the space variable cc £ C and the velocity 
V £ with d = 1, 2,3. Particles undergo the effect of collisions which are modelized 
here by a linear operator L acting on / through 

L{f){t,x,v) = p{t,x)Mfi{v) - fit,x,v), 

where 

p{t,x) = {f{t,x,v))( fit,x,v)dv. 


and where the equilibrium Mp is defined by 


Mp{v) = 


1 

(27r)'^/2 


exp {-\vf/2) 
m 



if /3 = 2 + d, 
if /3 £ (d, d + 2), 


( 1 . 1 ) 


where m is a normalization factor. In the sequel, we will always denote by brackets 
the integration over v. In order to capture a nontrivial asymptotic model, a suitable 
scaling has to be considered. The good scaling of the kinetic equation (see [24, 23, 2]) 
satisfied by the distribution function / is given by 


e°‘dtf + ev-\7a;f = L{f), with a =/? - d £ (0, 2], (1.2) 

where e: > 0 is the Knudsen number (ratio between the mean free path and a typical 
macroscopic length). The case a = 2 refers to the classical diffusion limit (the equi¬ 
librium corresponds to the first equation in (1.1)) whereas the case a £ (0,2) refers 
to the anomalous diffusion limit (the equilibrium corresponds to the second equation 
in (1.1)). Although our analysis could be applied to general operators L, we will 
consider for a sake of simplicity a BGK type operator. We suppose that Mp given by 

(1.1) is an even positive function such that {Mp{v)) = 1 and {vMp{v)) = 0. Equation 

(1.2) has to be supplemented with an initial condition /(O, x, v) = fo{x, v) and spatial 
periodic boundary conditions are considered. 

The main goal of this work is to construct numerical schemes for (1.2) which 
are uniformly stable along the transition from kinetic to macroscopic regime {i.e. for 
all e > 0). More precisely, two asymptotics are considered here: the diffusion limit 
(a = 2) and the anomalous diffusion limit (a £ (0, 2)). For the diffusion limit, several 
numerical schemes have already been proposed but in the anomalous diffusion limit, 
no numerical strategy has been proposed in the literature for (1.2), to the best of 
author’s knowledge. 

Mathematically, the derivation of diffusion type equation from kinetic equations 
such as (1.2) when a = 2 was first investigated in [29], [3], [21] or [12]. When Mp 
decreases quickly enough for large juj, the solution / of (1.2) converges, when e goes 
to zero, to an equilibrium function f{t,x,v) = p{t,x)Mp{v) where p(t,x) is solution 
of the diffusion equation 


dtp{t, a:) - Va: • {Dy^p{t, cc)) = 0, 


(1.3) 


and where the diffusion matrix D is given, in the simple case of the BGK operator, 
by the formula 


D= v®vMp{v)<lv. 


(1.4) 
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In particular, a crucial assumption on Mp for D to be finite is that 



\v\'^Mp(v)(iv < +oo; 


this is the case when Mp is Maxwellian (first case of (l-l))- 

Hence, if the equilibrium has no finite second order moment, the asymptotic of 
rescaled kinetic equation with a = 2 is not able to capture a nontrivial dynamics and 
the value of a should be tuned with Mp in order to recover this macroscopic equation 
when e 0. Typically, let Mp be a heavy-tailed distribution function, corresponding 
to the second case of (1.1). An example of such a distribution function we will often 
consider is Mp{v) = with m chosen such that {Mp{v)) = 1. 

In the case of astrophysical plasmas, it often occurs that the equilibrium Mp is a 
heavy-tailed distribution of particles which typically generates an anomalous diffusion 
behaviour (see [25, 28]). In particular, the diffusion matrix D is no longer well-defined 
and this requires to adapt the value of a in terms of the tail decreasing rate (3, in 
order to capture a nontrivial dynamics. It has been shown in [24, 23] and [2] that the 
appropriate scaling is a G (0,2) depends on the size /3 of the tail of the equilibrium 
through the relation fi = a + d. 

When e: goes to 0, the solution of this equation with BGK operator converges to 
p{t,x)Mp(v) where p is the solution of the anomalous diffusion equation which can 
be written in Fourier variable 


dtp{t,k) = -K\k\°‘p{t,k), 


(1.5) 


where p stands for the space Fourier variable of p, fc is the Fourier variable and k is 
a constant which depends only on the size of the tail of the equilibrium Mp{v) and is 
defined by 


/ 


{w ■ e)^ m 
R‘i 1 + (w • e)^ \w\^ 


dw. 


( 1 . 6 ) 


for any e G such that |e| = 1 (note that k does not depend on e). The anomalous 
diffusion equation can also be written in the space variable 


= (_A,)tp(i,a;), (1.7) 

where m is the normalization factor of Mp in (second case in l.I), F is the usual 
function defined by 

^- 1-00 

r(a;) = / t"-^e-*dt, 

Jo 

and Cd,a is a normalization constant given by 


^d,a. 


27r^+“F (I - f) ■ 


( 1 . 8 ) 


The fractional operator of the anomalous diffusion is easily defined by its Fourier 
transform 
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but has also an integral definition 


(-A^) = p{x) = Cd,aP-V- 


j 

Jr 




p{x + y) - p{x) 
--dy 


where P.V. denotes the principal value of the integral. 

The anomalous diffusion occurs not only in astrophysical plasmas but also in the 
study of granular media (see [14, 5, 4]), tokamaks (see [13]) or even in economy or 
social science: the Pareto distribution is a power tail distribution that satisfies the 
second case of (1.1). This kind of distribution is a common object used to modelise 
the repartition of sizes in a given set (asterods, cities, income, opinions...), see [20] 
for detailed explanations. At a microscopic level, these velocity distributions arise 
when the motion of the particles is no longer governed by a brownian process but by 
a Levy process (see [10], [11]). 

In this work, we want to construct a numerical scheme over the kinetic equation 
using a heavy-tailed distribution equilibrium, which is robust when the stiffness pa¬ 
rameter £ goes to 0. Asymptotic Preserving (AP) schemes already exist in the case 
of the diffusion limit where the equilibrium is a Maxwellian (see [16, 6, 7, 17, 19, 22, 
18, 27]) but the strategy used to obtain them cannot be rendered word for word to 
the anomalous diffusion case. As in the diffusion case, the difficulties when numer¬ 
ically solving the kinetic equation in the case of anomalous diffusion come from the 
stiff terms which require a severe condition on the numerical parameters. But in the 
anomalous diffusion case, there is an additional difficulty since one needs to take into 
account the large velocities arising in the equilibrium Mp. Indeed, the role of these 
high velocities is crucial and they have to be captured numerically in order to produce 
the anomalous diffusion operator when £ —>■ 0. 

To derive a numerical scheme which enjoys the AP property, different strategies 
are investigated in this work. The first one is based on a fully implicit scheme in 
time (both the transport and the collision term are considered implicit). Using this 
brute force strategy, the so-obtained numerical scheme is obviously AP in the diffusion 
limit but not in the anomalous diffusion one. Therefore, a suitable transformation of 
this implicit scheme is proposed in this work and leads to a scheme enjoying the AP 
property in the anomalous diffusion limit. 

From a computational point of view, a fully implicit scheme may be very expensive 
especially if one deals with high-dimensional problems. To avoid this implicitation, 
we propose another strategy which is based on a micro-macro decomposition. The 
distribution function / is written as a sum of an equilibrium part pMp plus a remain¬ 
der g. A model equivalent to (1.2) is then derived, composed of a kinetic equation 
satisfied by g and a macroscopic equation satisfied by p. The macroscopic flux is then 
reformulated to obtain a numerical scheme which enjoys the AP property and which 
does not require the implicitation of the transport term. 

The last strategy proposed in this work is based on a Duhamel formulation of 
the equation from which a AP numerical scheme is obtained; this approach bears 
similarities with [15]. This numerical scheme is efficient for the diffusion case and has 
to be adapted also to get an AP scheme in anomalous diffusion limit. Moreover, there 
are two main advantages in using such strategy. First we derive a closed equation 
on p and the computation of the full distribution function is not needed at any time. 
Second, this strategy ensures that the numerical scheme has uniform accuracy in 
£. Eventually, it would be easier to construct high order (in time) schemes on this 
formulation, this will be done in [9]. 
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The paper is organized as follows. In the next section we shall start with formal 
computations recalling how the diffusion and of anomalous diffusion limits can be 
obtained from the kinetic equation. Section 3 is devoted to the derivation of the 
different numerical schemes: a fully implicit scheme, a micro-macro scheme and a 
scheme based on a Duhamed formulation of (1.2). After the theoretical study of their 
properties, several numerical tests are presented in Section 4 in a one-dimensional case 
with periodic boundary conditions in space to illustrate and compare the efficiency of 
the various numerical schemes. 

2. A closed equation for the density and its diffusion asymptotics. This 
section is devoted to the derivation from (1.2) of a closed equation satisfied by p. 
This equation will be the basic equation from which we formally derive the diffusion 
and the anomalous diffusion equations. Of course, the formal derivation of diffusion 
and anomalous diffusion equations is classical but, for a sake of clarity, we present 
one way (among many others) to perform such derivation. In particular, the formal 
computation below will drive the construction of a class of our numerical schemes. 
We recall here the hypothesis on the equilibrium Mp{v) introduced in (1.1): it is an 
even function such that {Mp(v)) = 1, {vMp{v)) = 0. Then we will consider the two 
following cases mentioned in (1.1): a = 2 and 0 < a < 2 

Case 1. (v 0 vMp{v)) < -l-oo and a = 2. 

Case 2. {v ® vMp{v)) = oo and a G (0, 2). For example, Mp has a heavy tail of 
the form 

Tfl 

Mp{v) ^ /3 e (d, d-l-2), (2.1) 

|«|->oo \v\P 

with a positive normalization factor m ensuring {Mp) = 1. 

2.1. Space Fourier transform based computations. Starting from (1.2), 
we propose here a method based on spatial Fourier transform to formally derive the 
asymptotic equation in both diffusion and anomalous diffusion limits. It will also be 
the basis of the numerical methods we set in Section 3.3. 

Eventually, we remind that we denote by f{t, k, v) (resp. p{t, k)) 

f{t,k,v)= [ e~^’'""f{t,x,v)dx, 

the space Fourier transform of f{t,x,v) (resp. p{t,x)). We have 

Proposition 2.1. 

(i) Equation (1.2) formally implies the following exact equation on p 

t 

p{t, k) = k,v)^+ jj" e-^ Mp{v)) pit - e“s, k)ds. (2.2) 

(ii) In case 1 (diffusion scaling), when e —>• 0, p solves the diffusion equation 

dtpit, k) = -(^{k- vf Mp{v)^ pit, k). 

(hi) In case 2 (anomalous diffusion scaling), when e —>■ 0, p solves the anomalous 
diffusion equation 

dtpit,k) = -K\kff pit,k), 
with K defined by (1.6) and j5 — d = a. 
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Proof. 


(i) Formal derivation of the expression (2.2) from (1.2) 

We start by taking the Fourier transform of (1.2) 


dtf + ie^ °‘vkf = ^ i^pMp - /) . 


We solve in / (assuming p is given) 



which can be written after a change of variables s —>• (t — s)/£“ 



Still denoting by brackets the integration over the velocity space, we can integrate the 
previous expression with respect to v to get (2.2). Then, we can write it to make the 
limit equation appear 



—\esk-v 


Asp 


e 



(2.3) 


where 




(2.4) 


denotes the initial layer term. As we will see in the next parts of the proof, this 
formulation gives us the two asymptotic equations of the proposition. 

(ii) The diffusion limit 

Considering the diffusion scaling summarized in case 1, we will show that the limit 
equation is the diffusion equation. In the second integral of (2.3), we expand the 
expression of p in a power series of e up to the second order. Note that the first term 
of the expression is exponentially small. Indeed for all t > 0, Aq = o(e°^) when e goes 
to zero. Then we can write, using a Taylor expansion at £ = 0 for the exponential 


p = I e ^(MAdsd — ie I se ^(k-vMa)dsd -/ 


((fc • u)^Mp) dsp 



Reminding that {Mp) = 1 and that Mp is even, we can simplify this expression to 
obtain 


p = p-e^ {{k- v)'^Mp) p - e^dtp + o(£^), 


where we assume that p is sufficiently smooth in time to write 



Finally, when e goes to 0 we get the expected diffusion equation 


9tP = -{{k- v)‘^Mp) p. 
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(iii) The anomalous diffusion limit 

Now we will show that when Mp fulfills the conditions of case 2 (anomalous diffusion 
scaling), the solution of (1.2) converges to the solution of the anomalous diffusion 
equation when e goes to zero. 

As in the case of the diffusion limit, the first term has exponential decay when e goes 
to zero, Aq = o(e°°), for all t > 0. With the same assumption of smoothness in time 
for p, we get the following behaviour for the third term of (2.3) 


-a,p(t.k) + 0 ( 1 ). 


We now have to make appear the fractional diffusion in the second term of (2.3). 
Since the moment of order 2 of Mp is not finite, we cannot expand the exponential 
term into a power series in e anymore. We decompose this term as 


e-^ {e-'^^^-^Mp)ds= [^\-^ds+ e"" - l) ds, 

Jo Jo Jo 

and we focus on the second term on the right hand side of this last equality. 

In order to simplify the computations, we will consider the following specific form for 
the equilibrium Mp: Mp{v) = where m is chosen to ensure that {Mp(v)) = 

1. The general case, dealing with Mp satisfying (2.1) can be done similarly, see 
supplementary materials for details. We consider the integral 

^(e-iesfe'« _ 2 ^^ M/j) and, for nonzero k, we perform the change of variables w = e\k\v 


^-issk-v 


■» _ 1) _ 


m 


+ \v\P 


= (e 


\k\f-'^j (e-‘ 


-istct-ui 


TU'"' - 1 


{e\k\Y + \w\P 


-dw. 


Since the last integral has rotational symmetry, for any unitary vector e of we have 

= {e\k\Y~‘^C{s) + (e|A:|s)^“'^i?(£,s,/c), 


where 


C{s) = [ (e--- - 1) /^du;. 


( 2 . 6 ) 


and 


R{e,s,k) =—m{es\k\Y [ (e — l)- j- 

jR-i \w\P Y 


£:s|fc|)^ + \w\^ 


dw. 


In particular, R{e, s, k) tends to 0 when e go to zero and is bounded. If we set 
R{e, s, 0) = 0, equality (2.5) is also true for k = 0. 

From (2.3) and (2.5), we have to set a = /3 — d to get the fractional diffusion equation 
when e —> 0 and then we get 

p = p + e°‘\k\°‘ [ ( [ - l)e“®ds;0-e“9t;0 + o(e“), Vt > 0. 

Jo \JweR'‘ Pr / 
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From the equality 



1 ) 


;|/3 


dw 


Ms = - 


[ - 


(w ■ eY 


m 


gRd 1 + (w • e)2 \wY+< 


-dw := —K, 


we derive, when e goes to 0, the fractional diffusion equation (1.5) given in [24, 2]. 

□ 

2.2. Computations in the space variable. The aim of this section is to gen¬ 
eralize the previous computations to the original space variable. We have the following 
proposition 

Proposition 2.2. 

(i) Equation (1.2) formally implies the following exact equation on p 


p(t,x) = (e 



e (p{t — e°‘s,x — esv)Mp{v)) ds. 


(2.7) 


(ii) In case 1 (diffusion scaling), when e —>■ 0, p solves the diffusion equation 


dtp{t,x) = Va, • {DS/xPY^x)) , with D = {v® vMp(v)). 

(hi) In case 2 (anomalous diffusion scaling), when e —> 0, p solves the anomalous 
diffusion equation 


r, / , mr(a -f 1) , . , “ / X 

dtpft, x) = - - (-Aa;) 2 p(t, x), 

where F is the usual Euler Gamma function, m the normalization constant appearing 
in the expression of the equilibrium Mp, Cd,a is defined by (1.8) and (3 — d = a. 
Proof. 

(i) Formal derivation of the expression (2.7) from (1.2) 

As in the previous section, we can integrate (1.2) to get 

p = Aq-I- / {p[t — e°‘s,x — esv)Mp)ds, 

Jo 


with 


Ao{t,x) = (e J^foix-e^ °^tv,v)'j . (2.8) 

(ii) The diffusion limit 

Since (p(t, x — esv)) = {p(t, x + ssv)), we can rewrite (2.7) in order to make the limit 
equation appear 


p(t,x) = Ao{t,x) + / se 


p(t — £^s, X — esv) — p(t, X — esv) 


Mp{v) ) ds 


2 / Pit^x-evs) +p{t,x + evs)-2p{t,x) . 

Yf “ \ -JV- Mf{v}}ds 


+ e ^dsp{t,x). 


(2.9) 
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First, we again notice that the initial layer term Aq defined in (2.8) has exponential 
decay when e goes to zero. Second, Taylor expansions in the two integrals in time 
enable to write for e sufficiently small 


p = / se ® {-dtp) {Mp) ds + — s^e ® {^-vH^pvMp) ds 

+ fl-e ^ ' 


(l-e-^)p + o{e^), 


where H^p denotes the Hessian of p. This last expression can be simplified into 

dtp = {*vHxpvMp{v)) , 

when e goes to 0. With D defined in (1.4), we have {*vHxpvMp{v)) = • {DV^p) 

so this last equation is the diffusion equation as expected. 

(hi) The anomalous diffusion limit 
We rewrite (2.7) in the case of the anomalous diffusion scaling as follows 


p{t,x) = Ao{t,x) + £° 


.e- ) ds 


+ / e "" {{p{t,x - £vs) - p{t,x)) Mp{v)) ds + / e '‘ds {Mp{v)) p{t,x). 
Jo Jo 

( 2 . 10 ) 

Under smoothness assumptions on p, the first integral of (2.10) becomes, for small e 
_s / p{t — e'^SjX — evs) — p{t,x — evs) 


se 


£ S 


Mp{v) ) ds = -dtp{t, x) + o(l). 


Since the last term of (2.10) can be explicitly computed, let us focus on the second 
integral of (2.10). We perform the change of variables y = x — evs in the integral over 
V and exchanging the integrals in t and y to get 


Jo 

where 


((p(f, X — evs) — p{t, x)) AIp{v)) ds = 


/R<i 


p{i.y) - p{t,x) 
\x-yf 


a{e,x- y)dy, 


i{e,z) = 


(es)' 




ds. 


( 2 . 11 ) 


( 2 . 12 ) 


Now, we want to find an equivalent of a{e,z) for small values of e. Here we consider 
the equilibrium Mp{v) = the general case with Mp satisfying (2.1) can be done 

similarly (see supplementary materials for more details). The integral a becomes 


i{e,z) = ■ 


(es)^ 


) (es)^ + \z\P (es)^ 

and we consider, for a nonzero z, the following quantity 


ds, 


p+oo 

a{e, z) — me^“'^r(/3 — d + 1) = —m / (e. 

Jo 


sr ® 


(es)^ 


{es)P + \z\P 


ds 


^+oo 


— m 


{es)^ 


{es)P + |z|^ (es)' 


rds. 
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Then, for almost every z G (a(£,z) — me^ ‘^r(/3 — d+1)) converges to 0 

when £ goes to zero. We also have the following estimation 


i(e, z) - m£p-‘^T(P -d+l)\<K . 


(2.13) 


Let us remark that for a reasonable regularity of p, the integral appearing in (2.11) 
is well defined as a principal value because fd G (d, d + 2) 


P.V. 


/R'i 


p{t,y) - p{t,x) 
\y - xf 


dy = 


p{t, y) - p{t, x)-{y-x)- S/xpit, x) 




\y - x\P 


dy. 


Hence, using (2.13) we obtain that for almost all y G 


- <i + D) 

is dominated by an integrable function. We also obtained that for almost all y 
it tends to zero when £ goes to zero, by using the dominated convergence theorem for 
this function. 

We have, from (2.11) 


/o 


e ^ {{p{t,x — £vs) — p(t,x)) Mi 3 {v)) ds 

f p{t,y) - p{t,x) - {y - x) ■WxPit,x) 


V \x - y\^ 

where we used that a is even. So we have 


a{£,x- y)dy, 


f 


£—^0 


e^O 


e {{p{t,x — £vs) — p{t,x)) Mp{v)) ds 


m£^-‘^T{P-d+l) [ 

dR'i 

m£^~'^T{P — d + 1)P.14 j 
Jr'^ 


p{t, y) - p{t, x)-[y- x)Wxp{t, x) 


\y-x\f^ 

p{t,y) - p{t,x) 


dy 


\y - xf 


dy. 


We remark that we have to set /? — d = a to recover the limit equation, so we get 

pi.t,y) - pi.t,x) 


p{t,x) = —£°‘dtp{t,x)+m£°‘r{a + l)P.V. [ 
^^0 Jr 

that is when e goes to 0 

dtp{t,x) =—mT{a + 1)P.V. f 

Jr 


R<^ \y-x\<-+<i 

p{t,x) - p{t,y) 


dy + p(t,x) +o(e“). 


\y 


_ ^\a-\-d 


dy, 


which is the expected anomalous diffusion equation. 

□ 

3. Numerical schemes. This section is devoted to the presentation of appro¬ 
priate numerical schemes to capture the solution of (1.2). We want these schemes to 
be Asymptotic Preserving (AP), that is, for arbitrary initial condition /o 
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1. to be consistent with the kinetic equation (1.2) when the time step At tends 
to 0, with a fixed e, 

2. to degenerate into a scheme solving the asymptotic equation (diffusion or 
anomalous diffusion) when e goes to 0 with a fixed time step At. 

In the sequel, three different numerical schemes are presented: a fully implicit scheme, 
a micro-macro decomposition based scheme and a Duhamel based scheme. We will 
consider a time discretization t„ = nAt, n = 0,..., N such that NAt = T where T is 
the hnal time and we will set /" ~ /(tn). 

As we used Fourier transform in the previous formal computations, we will con¬ 
sider a bounded domain O for the spatial domain with periodic conditions. Hence, 
we will be able to use the discrete Fourier transform in the algorithm. 

3.1. Implicit scheme. The first idea to design a scheme for the kinetic equation 
(1.2) is to set an implicit scheme over the Fourier formulation of the kinetic equation. 
We will see that in the case of the diffusion limit it preserves easily the asymptotic 
equation. But in the case of the anomalous diffusion limit, the effect of large velocities 
must be taken into account to obtain the good asymptotic and this needs a suitable 
modification of this fully implicit scheme. 

We start with (1.2) written in the spatial Fourier variable and consider a fully 
implicit time discretization 


We have 


with 


- r 

At 


+ -(l + iefc-n)r+i = -p-+iM^. 


1 -I- lAefc ■ V 1-1- iXek ■ v 


A = 


At 


e°‘ + At' 


(3.1) 


(3.2) 


Note that 0 < A < 1, A —0 and that A —>■ 1. 

—^0 s —^0 

To compute from (3.1), has to be determined hrst. To do that we 
integrate (3.1) with respect to v to get 


that is 


;^n+l 


= A 


M, 


p 


1 -I- '\Xek ■ V 


/3"+i + (1 - A) 


r 


1 -I- iAefc • V 


P{TL-\-l 


1 


^1-1- iXek ■ V J 


Mfi 

l+iAefc-v 


1 / iXek-vM^ \ 

1 —A \ 1+iAefc-i; / 


(3.3) 


(3.4) 


We have the following proposition 

Proposition 3.1. In the case 1 (diffusion scaling), we consider the scheme 
defined for all k and for all time index 0 < n < N, NAt = T by 



Ma 


1 -I- iAefc • V . 


A 1 + iAefc • V 

(1 - A)r + 


/ iXek ■ vMp 


1 — A \ 1 -I- iXek ■ V 


(3.5) 



















12 


N. CROUSEILLES, H. HIVERT AND M. LEMOU 


with X = At/(e^ + At) and with the initial condition f^{k, v) = fo{k, v). This scheme 
has the following properties: 

(i) The scheme is of order 1 for any fixed e and preserves the total mass 

Vne Il,A],p'‘(0) = p°(0). 

(ii) The scheme is AP: for a fixed At, the scheme solves the diffusion equation 
when e goes to zero 


— jPfk) 
At 


= -{{k-vfMp{v))p^+^[k). 


(3.6) 


Proof. As (i) is straightforward, let us prove (ii). In the case of the diffusion 
asymptotic, we must check that the scheme is AP. We rewrite (3.4) to get 


(l-A)p”+‘ = - 


+ / 

\ 1 + iXek ■ V / \ 1 + 1 


r 


iXek ■ V 


that is, with the equality (t^) 


_ 1 
— 7^, 


pn+1 _ pn 

At 


leXk -v \ , 1 


1 + iAefc • V 


\Xek ■' 


At \ 1 + iXek ■ V 


r 


(3.7) 


Hence we can write 


^n+l _ 

At 


{k ■ v) X 


2 \2 


1 + X^£^{k ■ v)‘ 


:Mp ) - 


1 / iAefc • V 
At \ 1 + iAefc • v‘ 


and so, when e —>■ 0 with a fixed At, the scheme degenerates into 


^n+l _ -n 

At 


^=^-{{k-vfMp{v))p 


n+1 


that is an implicit scheme for (1.3) in Fourier variable. 

□ 

Now we consider the case of an anomalous diffusion scaling 0 < a < 2. In this 
case, the equilibrium Mp is heavy-tailed and D = {v ® vMp(v)) = -too. However, we 
have 


lim — 
e->o e“ 


(fc • V) A^ 

I + A^e^(fc • vf 


:Mpiv))p^+^ = \krK 


with K defined in (1.6). Unfortunately a discretization in velocity of the integral 
involved in this limit does not see the heavy tail of the equilibrium, which means that 
the effect of large velocities is completely missed. In fact, the limit when e goes to zero 
of the implicit kinetic scheme we wrote in Prop. 3.1 degenerates into the following 
scheme 
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which is not the correct asymptotics. Hence we have to transform its expression to 
make the right limit clearly appear in the scheme. The scheme is presented in the 
following proposition. 

Proposition 3.2. In the case 2 (anomalous diffusion scaling), we consider the 
following scheme defined for all k and for all time index 0 < n < N, NAt = T by 


r 


p-+\k) = 




1 +\\ek ■ V 


Mf, 


1 +\Xek ■ V 
1 


(gA|fc|)° 
1 - A 


{v ■ ef 


(eA 


1 + \\ek ■ V . 


(1 - A)r + 


+ \v\°^+d- 1 + (n- eY 


with A = At/(e“ + At), where e is any unitary vector and with the initial condition 
y(k,v) = fo{k,v). This scheme has the following properties: 

(i) The scheme is of order 1 for any fixed e and preserves the total mass 

Vne Il,iV],p"(0) = p°(0). 

(ii) The scheme is AP: for a fixed At, the scheme solves the anomalous diffusion 
equation when e goes to zero 


fp^^{k) — fP{k) 
At 


t^\kYfP+\k), 


(3.8) 


where k has been defined by (1-6). 

Proof As (i) is straightforward, let us prove (ii). We start by considering the 
expression of in the scheme of Prop. 3.1, and more precisely the integral 


ieAfc • V Mp 
1 + ieXk ■ V 


= / Mp{v) 


e^X? {k ■ vY 

-o di;. 

1 + e^A^ (fc • vY 


As in the previous section, we perform the change of variables w = eX\k\v to obtain 



er^A^ {k ■ vY 
1 + e^A^ (fc • vY 


dn = (eA|fc|) 




= (sAlfcl)--^ 




{w ■ eY 

1 + (w • eY 


dw, 


(3.9) 


with e = fc/|fc|; note that the last integral does not depend on this unitary vector e 
thanks to its rotational invariance. 

We consider explicitly the case Mp(y) = with m such that {Mp) = 1 and 

with a £]0, 2[. We have 


Mp 


e^X^ [k ■ v) 

1 + e^X^ {k ■ vY 


= (£A|fc|)^ 


m 


{w ■ eY 


K'i (eA|fc|)“''’‘^ + 1 + (w • ef 


rdw, 


where e is any unitary vector. Inserting this formula in (3.5), we get the expression for 
p"+i. The equation satisfied by in Prop. 3.1 is unchanged, then the numerical 
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scheme of the proposition is derived. Now, it remains to show that this scheme enjoys 
the AP property. From the following form of the equation on 

- P" ^ /_ - _ • e)^ \ .n+i _ / iAefc -v /" \ 

At \(eA|A:|)“+'' + |u|“+^l + (?;-e)7 \ 1 + iAeA: • u At / ’ 

(3.10) 

we easily observe that it degenerates when e goes to zero to an implicit scheme for 
the anomalous diffusion equation. This concludes the proof. 

□ 

3.2. Scheme based on a micro-macro decomposition. In the previous part, 
we constructed a fully implicit scheme enjoying the AP property. The implicit char¬ 
acter of the transport may induce a high computational cost. Therefore, we propose 
another scheme, which is based on a micro-macro decomposition of the kinetic equa¬ 
tion and in which the transport part is explicit in time. In the diffusion case, such a 
scheme has been set in [22]. We first recall here the way to derive it in the diffusion 
case. Then, we consider the case of the anomalous diffusion and show how to take 
into account the effect of the large velocities induced by the heavy-tailed structure of 
the equilibrium Mp{v). 

The distribution / is decomposed as 


f = if) Mp + g = pMp -h g, 


where (g) = 0. Injecting this decomposition into (1.2), we have 


dtpMf} + dtg + ■ VxpMp + -^v ■ V^g = -^9- (3.11) 

To derive an equation on p, we integrate with respect to v to obtain 

dtp+-^{vV,g) = 0. (3.12) 

Now we replace dtp in (3.11) to get an equation satisfied by g 


From this formulation the semi-implicit scheme proposed in [22] writes 

„n+l 


p 


9 


n+l 


At 


+ -{»-V,9"«> = 0 


+ At-. + A ( 9 . v,,," - {V . v,!,") Me) = -T, 


n+1 


(3.14) 


and we have the following proposition 

Proposition 3.3. In case 1 (diffusion scaling), we consider the following scheme 
defined for all x € fl, v € V and all time index 0 < n < N by 


0^+1 


1 


- + j{-V,9"«> = 0 

■ ^-P"Mp + -{v - {v ■ Mp) = 

Ai e e 


with initial conditions p^{x) = p(0,x) and g^(x,v) = f(0,x,v) — p(0, x)Mp(v). This 
scheme has the following properties: 
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(i) The scheme is of order 1 for any fixed e and preserves the total mass 

Vne Il,7V],p"(0) = p°(0). 

(ii) The scheme is AP: for a fixed At, the scheme solves the diffusion equation 
when e goes to zero 


p^+^{k) - p^{k) 

At 


= -{{k-vfMi3{v))p’^{k). 


For the proof, we refer to [22]. 

We now consider the anomalous diffusion scaling. Formulation (3.14) does not 
work because the fractional laplacian does not appear when e goes to 0, so we have 
to modify the micro-macro scheme. As Mp is even we have 

{v ■ V,g"+i) = {v • V^p'^+^Mp) + {v ■ = (u • V,r+i) , 

and so the first equation of (3.14) can be rewritten as 




At 


Now we use an implicit formulation of the kinetic equation to express 

_ jtn 


At 


+evv,r+^ = p^+^Mp-r+\ 


As we did in the previous part we introduce the variable A defined by (3.2) and write 
this last equation as 

(7 + eXv ■ V,) r+i = Xp^+^Mp + (1 - X)r, 


that is 

/"+! =X{I + sXv V,)-^ p^+^Mp + (1 - A) (7 + sXv ■ V,)’^ /". 

As A 0{At), we will use the following approximated expression for 

/-+! = A (7 + eAi; • V,)"^ p'^+^Mp + (1 - A) (7 - eXv ■ V,) /” + 0(At), 

to avoid a costly inversion of the transport operator. In this expression, the terms we 
removed also vanish when e goes to zero. The scheme presented in Prop. 3.3 can be 
rewritten in the anomalous diffusion limit as 

+ ^A (n • V, (7 + eXv ■ V,)”' p^+^Mp) 

+ ^(1 -X){v V, (7 - eXv ■ V,) D = 0, 
e“ 

+^v ■ V^p-^Mp + ^ (u • - {v ■ V, 5 ") Mp) = 


pU+l _ pn 

At 


gn+l _ gn 


At 


We can make additional simplifications. First, we remark that, 

{v • V, (7 - eXv ■ V,) n = {v ■ V,/") + 0{At) = {v ■ V,^") + 0{At) 


(3.15) 
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where we used {vMp) = 0 and A = 0{At) at fixed e. Then, the last term in the first 
line of the scheme (3.15) can be replaced by £^““(1 — X) {v ■ V Second, we will 
reformulate the second term of the first line of the scheme (3.15) (which produces the 
fractional diffusion when e goes to zero) in order to get the right limit when e 0. 
To simplify the presentation, we present the reformulation using the Fourier variable 
(because of the non-local character of the fractional laplacian in the space coordinates, 
the case in these coordinates is more delicate). The second term of (3.15) writes in 
the Fourier variable 



ik ■ V 

1 -I- ieXk ■ V 



p-+\k) 


1 


e^X^ {k ■ v)'^ 

1 -I- e^X^{k ■ vY 



p"+i(fc). 


Then, as in (3.9) we make the change of variables w = eA|fc|u in the velocity integration 
to get the following semi-discrete scheme (with Mp{v) = jq^^) 


p"+i - p" 

At 


e 


e“ -1- At 


|. |a / • e)' ^ 

' ' \l + (u-e)2(£A|fc|)/3 + |u|/3 
• V,p") = 0, 


5"+i(fc) 


^ • V.p”Mp + - (z;. - (z; • V,p") Mp) = 


(3.16) 


where e is any unitary vector and where we denoted by the inverse of the Fourier 
transform. The following proposition summarizes the main properties of this micro¬ 
macro scheme: 

Proposition 3.4. In case 2 (anomalous diffusion scaling), we introduce the 
following micro-macro the scheme defined for all x G v G and all time index 
0 < n < N, NAt = T by 


pU+l _ pn 

At 


e 


£“ -I- At 


|,|a / • e)' ^ 

\ 1 -I- (z; • e)^ {eX\k\)^ -|- |z)|^ 
{v ■ V,p") = 0 


;^n+l 


(k) 


_ ri^ pc- 1 

^ +-^ • V.p”Mp + - (z; • X7xg- - {v ■ X7x9n Mp) = --9^^\ 


where e is any unitary vector, and where initial conditions are p^{x) = p(0,a:) and 
g°{x,v) = f{0,x,v) — p{0,x)Mp. This scheme has the following properties: 

(i) The scheme is of order 1 for any fixed e and preserves the total mass 


Vne Il,A],p'‘(0) = p°(0). 


(ii) The scheme is AP: for a fixed At, the scheme solves the diffusion equation 
when e goes to zero 

where k is defined by (1-6). 

Proof The proof of (i) is immediate, let us prove (ii). We have when e goes 
to zero with a fixed At, the scheme solves the anomalous diffusion equation. The 
equation on g gives, when e goes to zero. 


p"+i = 0(min(£,£“)). 
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The equation on p gives, in the Fourier variable 


/5"+i - p" 


+ A“|fc|^ 


(v ■ e)^ 


At ''II \ 1-I-. e)2 _l_ |^|P 

Thanks to the following relations, 


;?n+l _ 


= C>(min(e^ei+“)). 


A“ = 


At 


£“ + At 


= 1 + 0(e°‘), and 


{eX\k\)f^ + \v\^ \v\P 


+ 0(£“+'^), 


we get the implicit discretization of the anomalous diffusion equation (3.8) when 
£ —>■ 0. Hence the scheme solves the anomalous diffusion equation when e goes to zero 
with a fixed At. □ 


3.3. Scheme based on an integral formulation of the equation. In the 

previous parts, we wrote two AP schemes solving (1.2) in the cases of the diffusion 
and anomalous diffusion limits, both of them were of order 1 in time. Here we present 
a scheme based on a Duhamel formulation of (1.2) which has uniform accuracy in 
£. This approach bears similarities with the UGKS scheme (see [30, 26]). As in the 
previous parts, the case of the classical diffusion gives easily an AP scheme but the 
large velocities effects require a specific treatment for the anomalous diffusion case. 
For both diffusion and anomalous diffusion, we start from (2.2) 

t 

pit, k) = k,v)^+ e-" (e-*"®'= "Mp(z;)) p(t - £“s, fc)ds, 

hence evaluating at time t = leads to 

p(trt+i; k) = Ao(fra-|-i, k) + / e (e ’'Mp) p(t„+i — £“s, fc)ds, (3.17) 

j=0 


where Aq is defined by (2.4). We perform a quadrature of order 2 in the integrals. 
Assuming that the time derivatives of p are uniformly bounded in £, we have: 

Vje [i,iV-i],Vse ■ 

/3(t„+i-£“s,/c) = a(£, s)p(t„+i-tj, fc) + (l-a(£, s))/3(t„+i-tj+i,fc) + 0(At^), (3.18) 

uniformly in £, with a(e,s) = 1 — Inserting (3.18) in the integral term of 

(3.17) leads to 

p{t„+i - £“s, k)ds = Cj(fc)/5(t„+i - tj,k) 

e“ 

+ bj{k)pitn+i - tj+i, k) + 0{AP) , 

(3.19) 


uniformly in £, where we used the following notations Vj C |0,A7] 

q+i 

bjik) = 


£“s — ti 


At 


Cjik) = 


U+i 


1 - 


£“s — tj 

aF 


e-‘""'="'A/p)ds, 

(3.20) 

M e-"(e-''^"'='’'Mp)ds. 

(3.21) 
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3.3.1. AP scheme for diffusion. We consider the case of the diffusion limit 
and we use the previous quadrature to write an AP scheme for the kinetic equation 
(1.2). We use the quadrature (3.19) in (3.17) to write 

n 

bjp(tn—j, fc) + O (At^)^ , (3.22) 

1=0 


and as nAt < T then ^ At^ = 0{At) and we get a first order scheme that writes 
1=0 


p^+\k) 


io(t„+i, fc) + ^ {c,p^+^-^{k) + b,p^-^{k)) + 6o/5"(fc) 
1=1 


1 - Co 


To ensure the AP property, the time integration in the last two terms are computed 
exactly to get 


b,{k) 


Cj{k) 


At (1 + eik ■ vY 


-(l+eifc-p) 


1 + e\k ■ 


Mp{v) 


g-'^(l+£ifc''u) 
1 + eik ■ V 


g—74r(l+eife-«) — g—^(l+eifc'«)^ 

At (1 + eik ■ vY 



(3.23) 

(3.24) 


We have the following proposition: 

Proposition 3.5. In case 1 (dijfusion scaling), we consider the following the 
scheme defined for all k and for all time index 0 < n < N such that NAt = T by 


io(t„+i, fc) + ^ {c,p^+^-Yk) + bjP^-Yk)) + bop^ik) 

Sn+l/7,', _ _ 1 = 1 _ 


p^+Yk) = 


1 - Co 


with the initial condition fP(k) = p(0,k) and Ao,bj and Cj defined in (2.4)-(3.23)- 
(3.24)- This scheme has the following properties: 

(i) The scheme is of order 1 and the total mass is conserved 

Vne Il,A],p'‘(0) = p°(0). 


(ii) The scheme is AP: for a fixed At, the scheme solves the diffusion equation 
when e goes to zero 

Remark 1. From the numerical test, the scheme appears to be of order one 
uniformly in epsilon, as highlighted in Fig. 4-'^- This uniform property can be proved; 
however, since the proof is rather technical, this is postponed to a work in preparation 
(see [9]). 

Proof From the construction of the scheme developed above, this scheme is clearly 
consistent and of order 1 in time, uniformly in e. By induction we also remark that 
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it conserves the numerical mass p"'(0). Let us prove that, for a fixed At the scheme 
solves the diffusion equation when e goes to zero. Indeed, for small e the terms 
bj(k),Cj{k) for j > 1 and Aq are exponentially small, hence the scheme becomes 


(1 — co{k))p'^'^^{k) = bo{k)p'^{k) + o(e°°), when e —>• 0. 


From the explicit expressions of 6o(fc) and co(fc) given in (3.23)-(3.24) we get 


bo{k) = ^ 

co(fc) = 1 - ^ • vfMp{v)) + o(s^). 

Hence when e goes to zero with a fixed At, the scheme of Prop. 3.5 degenerates into 


/5”+^(fc) — /3”(fc) 
At 


{{k-vfMp{v))p^+^{k), 


that is an implicit scheme solving the diffusion equation (1.3) in Fourier variable. 

□ 


3.3.2. AP scheme for anomalous diffusion. Now, we want to set an AP 
scheme for the kinetic equation (1.2) in the anomalous diffusion scaling. Still using 
the notations defined in (2.4)-(3.20)-(3.21), the scheme in Prop. 3.5 writes as 


(1 - co)/5"+i(fc) = io(Wi, fc) + E + bjp^-\k)) + 5o/3"(fc). 

i=i 


We remind that in the case of classical diffusion, the term co;o"+^ produces the right 
limit for small e. If we consider the same algorithm, as we are forced to consider 
a bounded velocity space, the numerical second order moment of Mp will be large 
but finite. Therefore, we have to suitably transform the integrals Cj and bj before 
discretizing in velocity. We have 


bj = 


£°‘s — t 


^e-*((e-‘^*'=''-l)M;3)ds+ I 


At 


e“s — tj 

-T —-e '^ds. 

At 


Ci = 


1 - 


e“s — L 


At 


I \ l)M/j)ds 
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and in the velocity integrations we set w = e\k\v as in (2.5), so when Mp{v) = 
m /(1 + kl^) we have 


h,= e‘^\kr 


+ 


m 


(e|fc|)^ + 


*j+i 

r 

^ T7T" 


- 1) ds 


e“s — tj 

At 


e ®ds, 


Cj= 


+ 


(e|fc|)^ + \w\^ 
e°‘s — t. 


1 - ) e"* - 1) ds 


1 - 


At 


1 1 


e ''ds, 


with e any unitary vector. To ensure the AP property, the time integrations are 
computed exactly to get 


hj=e 


*3 / 1 — e 


At 


_e-?- l-e“|fc| 




\ {e\k\)P + \v\^ 


(3.25) 


+e^\k\ 


m 


I 1 + if • e At (1 + iv ■ e)^ J {e\k\Y + \v\l^ / ’ 


1 - 


1 — e 


At 




m 


+£“|fc|^ 


0-'ji(l+i33-e) £« 


1 + iu • e At 




{e\k\)P + \v\P , 


(3.26) 


(1 + iu • e)^ j {e\k\Y + \v\P 


We have the following proposition 

Proposition 3.6. In case 2 (anomalous diffusion scaling), with the notations 
Aq, bj and cj defined in (2.4)-(3.25)-(3.26), we consider the following scheme defined 
for all k and for all time index 0 < n < N such that NAt = T by 


Ao{tn+uk) + ^ {cffp^+^-^{k) + bff^-^{k)) + 5o/5"(fc) 

^31+1/1,'i ___ 


p^+\k) = 


1 - Co 


with the initial condition pP(k) = p{0,k). This scheme has the following properties: 

(i) The scheme is of order 1 and preserves the total mass 

Vne Il,A],p''(0) = /i°(0). 

(ii) The scheme is AP: for a fixed At, the scheme solves the anomalous diffusion 
equation when £ goes to zero 


p”+^(/c) — fP(k) 

At 


= -n\k\^p^+ffk), 


where k has been defined by (1.6). 
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Remark 2. As in the case of diffusion, the numerical tests (see Fig. f.lO) 
suggests that this scheme is order one uniformly in epsilon. This uniform property 
can also he proved and this will he done in a future work (see [9]). 

Proof. As for the case of classical diffusion, the conservation of the mass is ob¬ 
tained by induction and the fact that it is of order 1 comes from the Taylor expansion 
we performed in the integrals. Let us prove that for a fixed At the scheme solves the 
anomalous diffusion equation when e goes to zero. From the exponential decay of Aq 
and bj,Cj,j > 1 we have when e goes to zero 

(l-co)/3"+i(fc)=6o(fc)/3"(fc)+o(e°°). 


On the one side we have 


bo{k) = ^ + o(e“), 


and on the other side we have 


co= 

= -£“|fc| 


1 


- 1 


1 + iv ■ e 
m {v ■ e)^ 


TTL \ 




= -e°‘\k\°‘K + l-— + oie°‘). 
' ' At ^ ' 


We deduce that when e goes to zero, the scheme degenerates into 


— /3”(fc) 
At 


-K\krp^+\k), 


where k is defined by (1.6). This is an implicit scheme for the anomalous diffusion 
equation (1.5). □ 

Remark 3. The previous schemes of Prop. 3.5-3.6 can he modified into a numeri¬ 
cal scheme which degenerates into a second order time approximation of the asymptotic 
model. Indeed, we can modify (3.22) to get the following scheme 


p'^+^k) = Ao(t„+i, fc) + 6o;5”(fc) + 5] + bjP^-Hk)) 

1=1 

+ + „„ - i)£!AA±£!W + ,1 _ 

By construction, this scheme degenerates when e goes to zero, into a Crank-Nicolson 
numerical scheme for the diffusion or anomalous diffusion equations. 

4. Numerical results. In this part, we present the numerical tests for the im¬ 
plicit scheme and the micro-macro scheme in the case of anomalous diffusion and for 
the scheme based on the Duhamel formulation of the kinetic equation in both cases 
of diffusion and of anomalous diffusion. In the sequel, we will denote by ISD (resp. 
ISA) the implicit schemes set in Prop. 3.1 (resp. in Prop. 3.2) in the diffusion case 
(resp. anomalous diffusion case); we will call MMSD (resp. MMSA) the micro-macro 
schemes in Prop. 3.3 (resp. in Prop. 3.4) in the diffusion case (resp. anomalous dif¬ 
fusion case). Eventually, the acronyms DSD and DSA will refer to the schemes based 
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on Duhamel Formulation in the Prop. 3.5 and 3.6. Finally, DS and ADS denote the 
numerical schemes (3.6) and (3.8) of the asymptotic models (diffusion and anomalous 
diffusion equations). 

In this part, we will consider t > 0,x € [—1,1],?; € R and the following initial 
data 


fo{x,v) = {l+sm{Trx))M/ 3 {v), 

and periodic boundary conditions are imposed in x. In the diffusion case, we will 
consider the equilibrium 


Mp{v) 



whereas in the anomalous diffusion case, we will take 


Mp{v) 


m 

1 + |?;|^ ’ 


with m chosen such that j^Mp{v)Av = 1. In the sequel, unless other discretizations 
are mentioned, the following numerical parameters will hold: 

(i) time: the final time will be set T = 0.1 and At = 10“^. 

(ii) space: we will consider a uniform mesh in x considering Nx = 64 points. 

(iii) velocity: in both diffusion and anomalous diffusion case, a uniform velocity 
grid of the truncated domain [—?;max, I'lnax] is considered, with Ny points. In the 
diffusion case, we will take Umax = 10 and Ny = 20. In the anomalous diffusion case, 
we will consider ?;max = 50 and Ny = 200. The discretization ensures at the discrete 
level that J^vMp{v)dv = 0. 

In the tests for the anomalous diffusion, we will consider a = 1.5. 

In the sequel, we will be interested in the relative error. The relative error is 
defined by the difference between a reference solution (obtained by a scheme using a 
refined mesh) and the solution obtained by the tested scheme through the following 
formula 


Error(T) 


11 Preference (T) Pscheme(T) |1 2 


IIPreference(T) |t2 


(4.1) 


where || • II 2 denotes the discrete norm and the reference scheme will be precised 
in each case. 


4.1. The implicit scheme in the case of the anomalous diffusion limit. 

In this section, we test the ISA scheme. We check that it is of order 1 for a fixed value 
of e and the AP property when e goes to 0 for a hxed value of At. 

Firstly, we remark that the ISD scheme used for the diffusion case does not 
preserve the asymptotic of anomalous diffusion: the left hand side of Fig. 4.1 shows 
that for £ = 10“® the implicit scheme returns a result very different from the result 
given by (3.8). Indeed, when £ « 1, the ISD scheme produces a stationary solution. 
So it appears that the change of variables in the velocity integrals is necessary to 
capture the anomalous diffusion limit, as shown in the right hand side of Fig. 4.1. 

The left hand side of Fig. 4.2 shows the error between the reference scheme (de¬ 
fined as the ADS scheme) and the ISA scheme as a function of £, for a = 0.8,1,1.5. 
We observe that the convergence of the kinetic equation to the anomalous diffusion 
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Fig. 4.1. For At = 10 the densities p{t = 0.1,x). Left: ISD seheme (with e = 10 by the 
ADS scheme and the initial data. Right: ISA scheme for different values of e and the ADS scheme. 


when e goes to zero arises with a speed a. Theorically, the convergence to the anoma¬ 
lous diffusion equation arises with a speed 2 -a)^ consider a finite 

domain for v in the tests, the convergence rate always appears to be e“. This result 
will be proved in [9]. 




Fig. 4.2. Left: For At = 10“^ and a = 0.8,1,1.5, the error as a function of e between the 
reference ADS scheme and the ISA scheme (log seale). Right: The error between the ISA scheme 
for a range of At = e°‘ and the ADS scheme as a function of At, which is fixed equal to e“ (log 
scale). 

However the convergence in time of the ISA scheme to the DSA scheme is not 
uniform with respect to e. Setting At = £“ in the scheme and letting At go to 0 shows 
that the densities does not converge to the density given by the anomalous diffusion 
equation. On the right hand side of Fig. 4.2, we plot the error in time between the 
ADS scheme as reference and the ISA scheme computed for a different time steps 
satisfying At = £“. We observe that the error does not converge to 0 when At goes 
to 0, illustrating the lack of uniformity. 

4.2. The micro-macro scheme. As we saw in Section 3.2, the MMSD scheme 
is known to be efficient in the diffusion case (see [22]). In this part, we focus on the 
MMSA scheme in the case of the kinetic equation with anomalous diffusion limit. 
As it does not use only the Fourier variable, we must consider a grid for the space 
variable. Here we will take x S [—1,1] discretized with = 32 points. As we 
solve the transport of the micro part (equation on g) in an explicit way (using finite 
differences), a CFL condition has to be imposed; then the time step is chosen small 
enough. 

We start by testing the consistency of the scheme. We fix a = 1.5 and for a 
decreasing range of time steps, we compute the solution given by the MMSA scheme. 
Then, for e = 1,0.5, we compare it to the reference solution given by MMSD with 
At = 10“®. For £ = 10“^, 10“®, we compare to the solution given by the ADS scheme. 
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Fig. 4.3 shows the error curves obtained for different values of e. It appears that the 
MMSA scheme is of order 1 in time. 




Fig. 4.3. For s = 1, 0.5,10 ^,10 the error in time between the MMSA scheme computed for 
a range of At and the reference scheme computed with At = 10~® (log scale). 


In order to show that the scheme preserves the asymptotic of anomalous diffusion 
we compute the densities p{t = 0.1, x) obtained by the MMSA scheme for a range of e; 
we compare them to the density obtained by the ADS scheme computed with the same 
At. In order to respect the CFL condition, we took At = 10“^ for e = 1, 0.5, 0.1, 0.01 
and At = 10“^ for the smallest e.The left hand side of Fig. 4.4 shows these densities 
in the case a = 1.5, we observe that they are very close to the anomalous diffusion 
limit for small e. On the right hand side of Fig. 4.4, the errors associated to this 
latter study are displayed in the cases a = 0.8, a = 1 and a = 1.5, the ADS scheme 
being considered as reference. We observe that the convergence happens with speed 
a, as expected. 




Fig. 4.4. The densities p{t = 0.1,3;). Left: MMSA scheme for different values of e and the DS 
scheme. Right: For At = 10“® and a = 0.8,1,1.5, the error with respect to e between the reference 
ADS scheme and the MMSA scheme (log scale). 


4.3. The integral formulation based scheme. In this section, we focus on the 
Duhamel formulation based schemes DSD and DSA in both diffusion and anomalous 
diffusion regimes. We put their properties of consistency, AP character and uniformity 
with respect to e, written in Prop. 3.5 and Prop. 3.6, in evidence in these two cases. 
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4.3.1. Case of diffusion limit. In the case of diffusion limit, we compare the 
results obtained by the DSD scheme to the results obtained with the ISD scheme. The 
first thing we check is the convergence of the algorithm when At goes to zero. We 
choose a range of e and we compare the densities obtained with the two algorithms 
as a function of At. Fig. 4.5 displays the errors associated to these comparisons. It 
appears that for e ^ 1, the scheme is of order 2 and of order I for small e. 




Fig. 4.5. For e = 1, 0.1,10 10 ®, the error in time between the DSD scheme and the reference 

ISD scheme computed for At = 10“® (log scale). 

In order to show the AP character of the scheme, we compute the densities given 
by the scheme DSD for a decreasing range of e (At being fixed) and we check that 
they converge to the solution of the diffusion equation as e goes to 0. The densities 
obtained are presented in Fig. 4.6. 



Fig. 4.6. For At = 10 ® the densities p(t = 0.1, a;) given by the DSD scheme for different 
values of e and the DS scheme. 

To highlight the fact that the scheme is first order in time uniformly in e, we 
compare the results given by the ISD scheme with At = I0“® to the results given by 
the DSD scheme for a range of At. These errors are displayed in the left hand side of 
Fig. 4.7 as a function of e where we observe that the error curves are stratihed with 
respect to At, showing the first order uniform accuracy of the scheme. This property 
appears also in the right hand side of Fig. 4.7 where we plotted the error between the 
DSD scheme computed for At = and the limit DS scheme computed for At = 10“^. 
This error tends to zero when At and e go to zero, thanks to the uniformity of the 
DSD scheme. 

4.3.2. Case of anomalous diffusion limit. In this section we test the DSA 
scheme. In the case of anomalous diffusion limit, we compare the results given by the 
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Fig. 4.7. Left: The difference between the ISD scheme computed for At = 10“® and the DSD 
scheme computed for a range of At as a function of s (log scale). Right: The difference between the 
DSD scheme computed for At = and the DS scheme computed for At = 10“^ as a function of 
At (log scale). 


DSA scheme to the results given by the DSA scheme computed with a smaller time 
step, to make the same properties appear. We start by highlighting the consistency 
of the scheme. Considering different values of e, we compute a reference solution with 
the scheme DSA with At = 10“® and we compare it to the results given by the DSA 
scheme for some larger time steps. Fig. 4.8 show error study of the convergence of 
the densities obtained. We observe that for e = 1,0.5 the scheme seems to be of order 
2 whereas it is of order 1 for smaller values of e. 



Fig. 4.8. Fore = 1,0.1,10 ^,10 the error as a function of At between the DSA scheme 
and the DSA scheme computed for At = 10“® (log scale). 


Then, we can study the AP character of the DSA scheme. The time step being 
fixed to At = 10“^, we check that the results given by the DSA scheme converge to 
the result given by the ADS scheme when e goes to 0. Fig. 4.9 presents the densities 
p{t = 0.1, x) obtained by these two schemes for different values of e and then, the 
error as a function of e is plotted for a = 0.8,1,1.5. As previously, the expected 
numerical speed of the convergence to the anomalous diffusion is recovered. 

To highlight the fact that the scheme is first order in time uniformly in e, we 
compare the results given by the DSA scheme with At = 10“^ to the results given 
by the same DSA scheme for a range of At. These errors are displayed in the left 
hand side os Fig. 4.10 as a function of e where we observe that the error curves are 
stratified with respect to At, showing the uniformity of the scheme with respect to e. 
As for the case of diffusion limit, the right hand side of Fig. 4.10 presents the error 
between the DSA scheme computed for At = e“ and the limit ADS scheme computed 
for At = 10“^. Since the DSA scheme is of order 1 uniformly in e, this error tends to 
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Fig. 4.9. Left: For At = 10“® the densities p{t = 0.1, x) given by the DSA scheme for different 
values of e and the ADS scheme. Right: For At = 10“® and a = 0.8,1,1.5, the error in e between 
the reference ADS scheme and the DSA scheme (log scale). 


zero when At and e go to zero. 




Fig. 4.10. Left: The difference between the DSA scheme computed for At = 10“® and the DSA 
scheme computed for a range of At as a function of s (log scale). Right: The difference between the 
DSA scheme computed for At = £°‘ and the DSA scheme computed for At = 10“^ as a function of 
At (log scale). 


5. Conclusion. In this paper, we have first presented a formal derivation of 
diffusion and anomalous diffusion asymptotics from a BGK kinetic equation. This 
analysis enables us to understand the role of the large velocities induced by the heavy¬ 
tailed equilibrium in the anomalous diffusion asymptotics. Moreover, this formal 
derivation paves the way to construct three different numerical schemes for the kinetic 
equation, which all enjoy the asymptotic preserving property in both the diffusion and 
anomalous diffusion regimes. 

The asymptotic preserving character in the case of anomalous diffusion comes 
from the fact that we managed to take into account the effect of large velocities of 
the equilibrium. As we saw in the simplest case of a direct fully implicit scheme, the 
asymptotic is not preserved when these large velocities are truncated. The schemes 
using the micro-macro formulation as well as the Duhamel formulation of the equa¬ 
tions require the use of the same trick on the velocities. Moreover, this last scheme 
seems to enjoy a uniform accuracy property which will be proved in [9]. 

In the near future, we aim at extending this work to more general context, con¬ 
sidering higher dimensions, non periodic boundary conditions. The case of singular 
collision frequency also may generate anomalous diffusion (see [1]) and this also de¬ 
serves a numerical study that we plan to do in a forthcoming work [8]. 
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